Análise dos dados de cancer colorretal - eventos recorrentes
Author
Camille Menezes
Published
July 3, 2024
Neste documento, duas abordagens, frequentista e Bayesiana, são consideradas para analisar o tempo até a ocorrência de múltiplas reinternações após a cirurgia de remoção do tumor de pacientes diagnosticados com câncer colorretal.
1 Análise descritiva
os dados de cancer colorretal apresentam 861 observações de 403. Desses indivíduos, 50.62% foram reinternados, 53.85% foram tratados com quimioterapia, 40.69% são do sexo feminino, 44.67% apresentam estágio tumoral A-B, 36.72% apresentam estágio tumoral C e 18.61% apresentam estágio tumoral D. É possível observar o Indice de Comorbidade de Charlson’s para esses indivíduos na figura a seguir.
Cerca de 80.14% das observações ocorreram até o terceiro evento. Para sintetizar os resultados obtidos, na análise descritiva foram consirados os tempos até a terceira internação.
Considerando todos os tempos, a primeira intenação ocorreu no 2º dia e a última no 2175º dia.
1.1 Tempo até o primeiro evento
Considerando o tempo até o primeiro evento, a primeira intenação ocorreu no 5º dia e a última no 2175º dia. As curvas de kaplan-Meier podem ser vistas no gráfico a seguir.
Valor-p do teste logrank
Valorp
chemo
0.026
sex
0.179
dukes
0.000
charlson
0.003
Pairwise comparisons using Log-Rank test
data: df1 and dukes
A-B C
C 0.013 -
D 1.2e-15 3.1e-06
P value adjustment method: bonferroni
Pairwise comparisons using Log-Rank test
data: df1 and charlson
0 1-2
1-2 1.0000 -
3 0.0015 1.0000
P value adjustment method: bonferroni
1.2 Tempo até o segundo evento
Considerando o tempo até o segundo evento, a primeira intenação ocorreu no 2º dia e a última no 1673º dia. As curvas de kaplan-Meier podem ser vistas no gráfico a seguir.
Valor-p do teste logrank
Valorp
chemo
0.196
sex
0.009
dukes
0.000
charlson
0.004
Pairwise comparisons using Log-Rank test
data: df2 and dukes
A-B C
C 0.45 -
D 1.5e-07 3.1e-05
P value adjustment method: bonferroni
Pairwise comparisons using Log-Rank test
data: df2 and charlson
0 1-2
1-2 0.021 -
3 0.043 1.000
P value adjustment method: bonferroni
1.3 Tempo até o terceiro evento
Considerando o tempo até o primeiro evento, a primeira intenação ocorreu no 8º dia e a última no 1572º dia. As curvas de kaplan-Meier podem ser vistas no gráfico a seguir.
Valor-p do teste logrank
Valorp
chemo
0.125
sex
0.093
dukes
0.000
charlson
0.018
Pairwise comparisons using Log-Rank test
data: df3 and dukes
A-B C
C 1 -
D 1.8e-08 1.3e-06
P value adjustment method: bonferroni
Pairwise comparisons using Log-Rank test
data: df3 and charlson
0 1-2
1-2 1.000 -
3 0.013 0.903
P value adjustment method: bonferroni
2 Modelo completo
2.1 Abordagem frequentista
A partir da abordagem frequentista, é possível notar que todas as covariáveis, com exceção da quimioterapia, são significativas ao nível de 5%. É possível notar ainda que, ao nível de 5%, há evidências o suficiente para concluir que a variância da fragilidade não é nula e o modelo não se reduz ao modelo para dados independentes.
coef
exp(coef)
exp(-coef)
se(coef)
z
p-value
2.5%
97.5%
Intercepto
-6.888
0.001
980.028
0.306
-22.515
0.000
-7.487
-6.288
chemoTreated
-0.012
0.988
1.012
0.119
-0.104
0.917
-0.246
0.221
sexFemale
-0.380
0.684
1.462
0.114
-3.320
0.001
-0.604
-0.156
dukesC
0.311
1.364
0.733
0.134
2.319
0.020
0.048
0.573
dukesD
1.150
3.159
0.317
0.167
6.883
0.000
0.823
1.478
charlson1-2
0.606
1.832
0.546
0.220
2.757
0.006
0.175
1.036
charlson3
0.177
1.193
0.838
0.124
1.426
0.154
-0.066
0.420
gama
0.906
2.475
0.404
0.039
23.207
0.000
0.830
0.983
xi
0.145
1.156
0.865
0.070
2.082
0.037
0.009
0.282
Tomando o exponencial do seu coeficiente estimado, mantidas fixas as outras covariáveis, tem-se que:
O risco de reinternação dos pacientes do sexo feminino é 0.684 vezes o risco dos pacientes do sexo masculino;
O risco de reinternação dos pacientes os pacientes que apresentam estágio tumoral de Dukes A-B é 0.733 vezes o risco dos pacientes que apresentam estágio C;
O risco de reinternação dos pacientes os pacientes que apresentam estágio tumoral de Dukes A-B é 0.317 vezes o risco dos pacientes que apresentam estágio D;
O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é 0.546 vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 1-2;
O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é 0.838 vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 3.
2.2 Resíduos
2.3 Abordagem Bayesiana
Utilizando o algoritmo Adaptative Metropolis-within-Gibbs, obteve-se uma taxa de aceitação de 43.95%, próxima do considerado ótimo.
Tomando o exponencial do seu coeficiente estimado, mantidas fixas as outras covariáveis, tem-se que:
O risco de reinternação dos pacientes do sexo feminino é 0.686 vezes o risco dos pacientes do sexo masculino;
O risco de reinternação dos pacientes os pacientes que apresentam estágio tumoral de Dukes A-B é 0.76 vezes o risco dos pacientes que apresentam estágio C;
O risco de reinternação dos pacientes os pacientes que apresentam estágio tumoral de Dukes A-B é 0.311 vezes o risco dos pacientes que apresentam estágio D;
O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é 0.526 vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 1-2;
O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é 0.827 vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 3;
Tabela: Medidas resumo das distribuições dos parâmetros da distribuição a posteriori, do desvio e do LP.
mode
exp.coef.
exp..coef.
Mean
SD
MCSE
ESS
LB
Median
UB
beta[1]
-6.949
0.001
1042.107
-6.944
0.317
0.014
924.227
-7.522
-6.944
-6.333
beta[2]
-0.015
0.985
1.015
-0.018
0.124
0.003
2000.000
-0.259
-0.018
0.225
beta[3]
-0.377
0.686
1.458
-0.385
0.120
0.003
2000.000
-0.621
-0.385
-0.152
beta[4]
0.275
1.317
0.76
0.302
0.138
0.003
2000.000
0.032
0.299
0.581
beta[5]
1.167
3.212
0.311
1.161
0.172
0.004
1792.647
0.820
1.165
1.491
beta[6]
0.643
1.902
0.526
0.607
0.230
0.005
2000.000
0.139
0.618
1.023
beta[7]
0.19
1.209
0.827
0.179
0.126
0.003
2000.000
-0.072
0.182
0.421
gama
0.915
2.497
0.401
0.914
0.041
0.002
904.900
0.836
0.916
0.993
xi
0.169
1.184
0.845
0.191
0.080
0.002
1615.365
0.060
0.183
0.372
Deviance
7403.471
-
-
7404.888
4.723
0.131
1622.999
7398.411
7404.078
7415.384
LP
-
-
-
-3737.795
2.363
0.065
1624.695
-3743.060
-3737.397
-3734.552
A partir dos gráficos a seguir e da Tabela, é possível notar que as distribuições a posteriori dos betas, gama e xi são simétricas. É possível notar também há uma notável estacionariedade, evidenciada pela ausência de tendências nos gráficos. O que indica que as amostras estão sendo colhidas de forma aleatória da mesma distribuição alvo.
Os correlogramas exibidos revelam uma escassa presença de autocorrelações significativas, sugerindo que as amostras podem ser consideradas independentes. Ou seja, as amostras geradas não são fortemente influenciadas pelas amostras geradas anteriormente.
Os desvios sugerem que a maioria das verossimilhanças foram maximizadas durante o processo, já que a distribuição é assimétrica à direita. Além disso, o correlograma não apresenta indícios de autocorrelação nos desvios das amostras.
Os logaritmos da distribuição a posteriori, das amostras indicam a sua boa adequação, uma vez que a sua distribuição é assimétrica à esquerda e o correlograma apresenta apresenta apenas três correlações significativas.
A partir da abordagem frequentista, é possível notar que todas as covariáveis são significativas ao nível de 5%. É possível notar ainda que, ao nível de 5%, há evidências o suficiente para concluir que a variância da fragilidade não é nula e o modelo não se reduz ao modelo para dados independentes.
coef
exp(coef)
exp(-coef)
se(coef)
z
p-value
2.5%
97.5%
Intercepto
-6.895
0.001
987.494
0.290
-23.767
0.000
-7.464
-6.327
sexFemale
-0.379
0.684
1.461
0.114
-3.318
0.001
-0.603
-0.155
dukesC
0.315
1.371
0.730
0.126
2.504
0.012
0.069
0.562
dukesD
1.155
3.173
0.315
0.162
7.116
0.000
0.837
1.473
charlson1-2
0.602
1.826
0.548
0.219
2.754
0.006
0.174
1.031
charlson3
0.176
1.193
0.838
0.124
1.426
0.154
-0.066
0.419
gama
0.906
2.474
0.404
0.039
23.258
0.000
0.830
0.982
xi
0.145
1.156
0.865
0.070
2.082
0.037
0.009
0.281
Tomando o exponencial do seu coeficiente estimado, mantidas fixas as outras covariáveis, tem-se que:
O risco de reinternação dos pacientes do sexo feminino é 0.684 vezes o risco dos pacientes do sexo masculino;
O risco de reinternação dos pacientes os pacientes que apresentam estágio tumoral de Dukes A-B é 0.73 vezes o risco dos pacientes que apresentam estágio C;
O risco de reinternação dos pacientes os pacientes que apresentam estágio tumoral de Dukes A-B é 0.315 vezes o risco dos pacientes que apresentam estágio D;
O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é 0.548 vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 1-2;
O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é 0.838 vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 3.
3.2 Resíduos
3.3 Abordagem Bayesiana
Utilizando o algoritmo Adaptative Metropolis-within-Gibbs, obteve-se uma taxa de aceitação de 44.07%, próxima do considerado ótimo.
Tomando o exponencial do seu coeficiente estimado, mantidas fixas as outras covariáveis, tem-se que:
O risco de reinternação dos pacientes do sexo feminino é 0.693 vezes o risco dos pacientes do sexo masculino;
O risco de reinternação dos pacientes os pacientes que apresentam estágio tumoral de Dukes A-B é 0.72 vezes o risco dos pacientes que apresentam estágio C;
O risco de reinternação dos pacientes os pacientes que apresentam estágio tumoral de Dukes A-B é 0.31 vezes o risco dos pacientes que apresentam estágio D;
O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é 0.536 vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 1-2;
O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é 0.837 vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 3;
Tabela: Medidas resumo das distribuições dos parâmetros da distribuição a posteriori, do desvio e do LP.
mode
exp.coef.
exp..coef.
Mean
SD
MCSE
ESS
LB
Median
UB
beta[1]
-6.899
0.001
991.283
-6.968
0.295
0.011
1093.123
-7.573
-6.966
-6.403
beta[2]
-0.367
0.693
1.443
-0.385
0.118
0.003
2000.000
-0.619
-0.379
-0.155
beta[3]
0.328
1.388
0.72
0.317
0.131
0.003
2000.000
0.060
0.317
0.574
beta[4]
1.172
3.228
0.31
1.167
0.170
0.005
1709.159
0.846
1.165
1.506
beta[5]
0.624
1.866
0.536
0.597
0.223
0.005
2000.000
0.135
0.606
1.007
beta[6]
0.178
1.195
0.837
0.180
0.124
0.003
2000.000
-0.067
0.181
0.422
gama
0.915
2.497
0.401
0.916
0.040
0.001
1028.286
0.839
0.915
0.996
xi
0.162
1.176
0.85
0.189
0.076
0.002
1803.806
0.067
0.182
0.353
Deviance
7402.034
-
-
7403.655
4.033
0.097
2000.000
7397.782
7403.005
7412.654
LP
-
-
-
-3732.806
2.019
0.048
2000.000
-3737.296
-3732.489
-3729.860
A partir dos gráficos a seguir e da Tabela, é possível notar que as distribuições a posteriori dos betas, gama e xi são simétricas. É possível notar também há uma notável estacionariedade, evidenciada pela ausência de tendências nos gráficos. O que indica que as amostras estão sendo colhidas de forma aleatória da mesma distribuição alvo.
Os correlogramas exibidos revelam uma escassa presença de autocorrelações significativas, sugerindo que as amostras podem ser consideradas independentes. Ou seja, as amostras geradas não são fortemente influenciadas pelas amostras geradas anteriormente.
Os desvios sugerem que a maioria das verossimilhanças foram maximizadas durante o processo, já que a distribuição é assimétrica à direita. Além disso, o correlograma não apresenta indícios de autocorrelação nos desvios das amostras.
Os logaritmos da distribuição a posteriori, das amostras indicam a sua boa adequação, uma vez que a sua distribuição é assimétrica à esquerda e o correlograma apresenta apresenta apenas três correlações significativas.
Na abordagem frequentista, o teste da razão de verossimilhanças resultou num valor-p de 0.9189116. Dessa forma, não rejeitamos a hipótese nula ao nível de 5%. Ou seja, há evidências o suficiente para concluir que o modelo reduzido é mais adequado que o modelo completo. Além disso, o AIC do modelo reduzido é menor, indicando a sua melhor adequação.
Tabela: AIC e BIC dos modelos completo e reduzido para os dados de cancer colorretal.
Modelo completo
Modelo reduzido
AIC
7409.446
7407.456
BIC
7429.237
7429.247
Na abordagem Bayesiana, o DIC para o modelo reduzido é menor que o do modelo completo, indicando a sua melhor adequação.
Tabela: DIC das amostras a posteriori dos modelos completo e reduzido para os dados de cancer colorretal.
Modelo completo
Modelo reduzido
Dbar
7404.89
7403.66
pd
11.15
8.13
DIC
7416.04
7411.79
O fator de Bayes é 0.88, isso indica que há evidências a favor do modelo reduzido. Além disso, o modelo com todas as covariáveis precisou de 17.02h para gerar todas as amostras, enquanto o modelo sem a covariável quimioterapia precisou de 14.55h.
Source Code
---title: "Análise dos dados de cancer colorretal - eventos recorrentes"author: "Camille Menezes"date: "2024/07/03"editor: sourceformat: html: number-sections: true smooth-scroll: true link-external-newwindow: true max-width: 3000px toc: true toc-location: left self-contained: true toc-title: "Sumário" page-layout: article theme: light: cosmo dark: darkly code-fold: show code-tools: true code-block-border-left: falses code-line-numbers: false code-copy: true code-overflow: wrapexecute: warning: false message: false echo: false---```{r}set.seed(25102022)#set.seed(20231120)``````{r}library(modeest)library(dplyr)library(gridExtra)library(ggpubr)library(tidyverse)library(survival)library(flexsurv)library(plotly)library(zoo)library(patchwork)library(knitr)library(forecast)library(HelpersMG)library(survminer)library(qqplotr)library(EstimationTools)library(LaplacesDemon)``````{r}df <-read.csv("readmission.csv")|>mutate(chemo =relevel(factor(chemo), "NonTreated"),sex =relevel(factor(sex), "Male"),dukes =relevel(factor(dukes), "A-B"),charlson =relevel(factor(charlson), "0"),death =NULL,time = t.stop) ```Neste documento, duas abordagens, frequentista e Bayesiana, são consideradas para analisar o tempo até a ocorrência de múltiplas reinternações após a cirurgia de remoção do tumor de pacientes diagnosticados com câncer colorretal.# Análise descritiva```{r}df1 <- df |>filter(enum ==1)df2 <- df |>filter(enum ==1& event ==0| enum ==2& event ==0| enum ==2& event ==1)df3 <- df |>filter(enum ==1& event ==0| enum ==2& event ==0| enum ==3& event ==0| enum ==3& event ==1)df4 <- df |>filter(enum ==1& event ==0| enum ==2& event ==0| enum ==3& event ==0| enum ==4& event ==0| enum ==4& event ==1)df5 <- df |>filter(enum ==1& event ==0| enum ==2& event ==0| enum ==3& event ==0| enum ==4& event ==0| enum ==5& event ==0| enum ==5& event ==1)df6 <- df |>filter(enum ==1& event ==0| enum ==2& event ==0| enum ==3& event ==0| enum ==4& event ==0| enum ==5& event ==0| enum ==6& event ==0| enum ==6& event ==1)```os dados de cancer colorretal apresentam `r nrow(df)` observações de `r length(unique(df$id))`. Desses indivíduos, `r round(1-table(summarise(group_by(df,id),sum(event))[,2])[1]/length(unique(df$id)),4)*100`% foram reinternados,`r as.numeric(round(prop.table(table(summarise(group_by(df,id),unique(chemo)) [2]))[2]*100,2))`\% foram tratados com quimioterapia, `r as.numeric(round(prop.table(table(summarise(group_by(df,id),unique(sex)) [2]))[2]*100,2))`\% são do sexo feminino, `r round(prop.table(table(summarise(group_by(df,id),unique(dukes))[2]))*100,2)[1]`% apresentam estágio tumoral A-B, `r round(prop.table(table(summarise(group_by(df,id),unique(dukes))[2]))*100,2)[2]`% apresentam estágio tumoral C e `r round(prop.table(table(summarise(group_by(df,id),unique(dukes))[2]))*100,2)[3]`% apresentam estágio tumoral D. É possível observar o Indice de Comorbidade de Charlson's para esses indivíduos na figura a seguir.```{r}d2 <- df6 |>group_by(enum,charlson) |>summarize(count =n()) |>mutate(pct = count/sum(count)) ggplot(d2, aes(fill=charlson, y=pct,x=factor(enum))) +geom_col(position =position_dodge(0.7),width =0.5) +scale_y_continuous(labels=scales::percent, breaks=seq(0,1,0.2))+scale_fill_manual(values=c("#5639cc","#26a4d0","#071e26"))+labs(y="Porcentagem de pacientes", fill="Índ. de Comorbidade",x="Evento")+theme_bw() ```Cerca de `r round(sum(prop.table(table(df$enum))[1:3])*100,2)`% das observações ocorreram até o terceiro evento. Para sintetizar os resultados obtidos, na análise descritiva foram consirados os tempos até a terceira internação.::: panel-tabset### TTT plot```{r}plot(TTTE_Analytical(Surv(t.start,t.stop, event)~1,method="cens",data = df))```### Pacientes versus tempos```{r}n <-c(1:43)d <-summarise(group_by(df[n,],id),max(time),) |>mutate(id=as.character(id)) |>data.frame()ggplot(df[n,])+aes(x=time, y=id, color=factor(event))+geom_point()+geom_segment(y=d[1,1],x=0,xend=d[1,2],color="#071e26")+geom_segment(y=d[2,1],x=0,xend=d[2,2],color="#071e26")+geom_segment(y=d[3,1],x=0,xend=d[3,2],color="#071e26")+geom_segment(y=d[4,1],x=0,xend=d[4,2],color="#071e26")+geom_segment(y=d[5,1],x=0,xend=d[5,2],color="#071e26")+geom_segment(y=d[6,1],x=0,xend=d[6,2],color="#071e26")+geom_segment(y=d[7,1],x=0,xend=d[7,2],color="#071e26")+geom_segment(y=d[8,1],x=0,xend=d[8,2],color="#071e26")+geom_segment(y=d[9,1],x=0,xend=d[9,2],color="#071e26")+geom_segment(y=d[10,1],x=0,xend=d[10,2],color="#071e26")+geom_segment(y=d[11,1],x=0,xend=d[11,2],color="#071e26")+geom_segment(y=d[12,1],x=0,xend=d[12,2],color="#071e26")+geom_segment(y=d[13,1],x=0,xend=d[13,2],color="#071e26")+geom_segment(y=d[14,1],x=0,xend=d[14,2],color="#071e26")+geom_segment(y=d[15,1],x=0,xend=d[15,2],color="#071e26")+geom_segment(y=d[16,1],x=0,xend=d[16,2],color="#071e26")+geom_segment(y=d[17,1],x=0,xend=d[17,2],color="#071e26")+geom_segment(y=d[18,1],x=0,xend=d[18,2],color="#071e26")+geom_segment(y=d[19,1],x=0,xend=d[19,2],color="#071e26")+geom_segment(y=d[20,1],x=0,xend=d[20,2],color="#071e26")+labs(y="Pacientes",x="Tempos até as reinternações ou censura (dias)",color="Evento")+scale_y_discrete(limits=as.character(1:20))+scale_color_manual(values =c("0"="#26a4d0","1"="#071e26"),labels =c("Censura", "Falha"))+theme_bw()```### Todos os pacientes vesus tempos```{r}ggplot(df)+aes(x=time, y=id, colour=factor(event))+geom_point()+labs(y="Pacientes",x="Tempos",color="Evento")+scale_colour_manual(values =c("0"="#26a4d0","1"="#071e26"),labels =c("Censura", "Falha"))+theme_bw()```### Índ. versus tempos```{r}`Tempo até a primeira reinternação`<- df$time Paciente <-index(df$time) Status <-ifelse(as.character(df$event) =="0", "Censura", "Falha") d3 <-data.frame(`Tempo até a primeira reinternação`, Paciente, Status)ggplotly(ggplot(d3)+aes(x=`Tempo até a primeira reinternação`, y=Paciente, colour = Status)+geom_point()+scale_colour_manual(values =c("Censura"="#26a4d0","Falha"="#071e26"),labels =c("Censura", "Falha"))+labs(x="Tempos", y="Índice", color="")+theme_bw())```::: Considerando todos os tempos, a primeira intenação ocorreu no `r min(subset(df, df$event == 1)$time)`º dia e a última no `r max(subset(df, df$event == 1)$time)`º dia.## Tempo até o primeiro eventoConsiderando o tempo até o primeiro evento, a primeira intenação ocorreu no `r min(subset(df1, df1$event == 1)$time)`º dia e a última no `r max(subset(df1, df1$event == 1)$time)`º dia. As curvas de kaplan-Meier podem ser vistas no gráfico a seguir.```{r}ekm1 <-survfit(Surv(time,event)~chemo, data=df1)tc1 <- ekm1$time[1:161]tc2 <- ekm1$time[162:(length(ekm1$time)-1)]sc1 <- ekm1$surv[1:161]sc2 <- ekm1$surv[162:(length(ekm1$time)-1)]d1 <-data.frame(Tempo =c(tc1,tc2), `Sobrevivência estimada`=round(c(sc1,sc2),2), Quimioterapia=c(rep("Não tratado",length(tc1)),rep("Tratado",length(tc2))))fig1 <-ggplot(d1)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Quimioterapia))+scale_color_manual(values=c("#071e26","#26a4d0"))+labs(x="", y="")+theme_bw()ekm2 <-survfit(Surv(time,event)~sex, data=df1)ts1 <- ekm2$time[1:192]ts2 <- ekm2$time[193:(length(ekm2$surv)-1)]ss1 <- ekm2$surv[1:192]ss2 <- ekm2$surv[193:(length(ekm2$surv)-1)]d2 <-data.frame(Tempo =c(ts1,ts2), `Sobrevivência estimada`=round(c(ss1,ss2),2), Sexo=c(rep("Masculino",length(ts1)),rep("Feminino",length(ts2))))fig2 <-ggplot(d2)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Sexo))+scale_color_manual(values=c("#071e26","#39cccc"))+labs(x="", y="")+theme_bw()ekm3 <-survfit(Surv(time,event)~dukes, data=df1)td1 <- ekm3$time[1:152]td2 <- ekm3$time[154:289]td3 <- ekm3$time[290:length(ekm3$time)]sd1 <- ekm3$surv[1:152]sd2 <- ekm3$surv[154:289]sd3 <- ekm3$surv[290:length(ekm3$surv)]d3 <-data.frame(Tempo =c(td1,td2,td3), `Sobrevivência estimada`=round(c(sd1,sd2,sd3),2), Dukes=c(rep("A-B",length(td1)),rep("C",length(td2)),rep("D", length(td3))))fig3 <-ggplot(d3)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Dukes))+scale_color_manual(values=c("#5639cc","#071e26","#26a4d0"))+labs(x="", y="")+theme_bw()ekm4 <-survfit(Surv(time,event)~charlson, data=df1)tcc1 <- ekm4$time[1:235]tcc2 <- ekm4$time[237:252]tcc3 <- ekm4$time[253:length(ekm4$time)]scc1 <- ekm4$surv[1:235]scc2 <- ekm4$surv[237:252]scc3 <- ekm4$surv[253:length(ekm4$surv)]d4 <-data.frame(Tempo=c(tcc1,tcc2,tcc3), `Sobrevivência estimada`=round(c(scc1,scc2,scc3),2), Charlson =c(rep("0",length(tcc1)),rep("1-2",length(tcc2)),rep("3", length(tcc3))))fig4 <-ggplot(d4)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Charlson))+scale_color_manual(values=c("#5639cc","#071e26","#26a4d0"))+labs(x="", y="")+theme_bw()``````{r}ggarrange(fig1,fig2,fig3,fig4) |>annotate_figure(bottom=text_grob("Tempo até a primeira reinternação (dias)"),left=text_grob("S(t)"))``````{r}valorp <-c(survdiff(Surv(time,event)~chemo, data=df1)$pvalue,survdiff(Surv(time,event)~sex, data=df1)$pvalue,survdiff(Surv(time,event)~dukes, data= df1)$pvalue,survdiff(Surv(time,event)~charlson, data= df1)$pvalue)d <-data.frame(Valorp =round(valorp,3))row.names(d) <-c("chemo", "sex", "dukes", "charlson")kable(d,caption ="Valor-p do teste logrank")``````{r}pairwise_survdiff(Surv(time,event)~dukes, data= df1, p.adjust.method ="bonferroni")pairwise_survdiff(Surv(time,event)~charlson, data= df1, p.adjust.method ="bonferroni")```## Tempo até o segundo eventoConsiderando o tempo até o segundo evento, a primeira intenação ocorreu no `r min(subset(df2, df2$event == 1)$time)`º dia e a última no `r max(subset(df2, df2$event == 1)$time)`º dia. As curvas de kaplan-Meier podem ser vistas no gráfico a seguir.```{r}ekm1 <-survfit(Surv(time,event)~chemo, data=df2)tc1 <- ekm1$time[1:173]tc2 <- ekm1$time[174:(length(ekm1$time)-1)]sc1 <- ekm1$surv[1:173]sc2 <- ekm1$surv[174:(length(ekm1$time)-1)]d1 <-data.frame(Tempo =c(tc1,tc2), `Sobrevivência estimada`=round(c(sc1,sc2),2), Quimioterapia=c(rep("Não tratado",length(tc1)),rep("Tratado",length(tc2))))fig1 <-ggplot(d1)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Quimioterapia))+scale_color_manual(values=c("#071e26","#26a4d0"))+labs(x="", y="")+theme_bw()ekm2 <-survfit(Surv(time,event)~sex, data=df2)ts1 <- ekm2$time[1:221]ts2 <- ekm2$time[222:(length(ekm2$surv)-1)]ss1 <- ekm2$surv[1:221]ss2 <- ekm2$surv[222:(length(ekm2$surv)-1)]d2 <-data.frame(Tempo =c(ts1,ts2), `Sobrevivência estimada`=round(c(ss1,ss2),2), Sexo=c(rep("Masculino",length(ts1)),rep("Feminino",length(ts2))))fig2 <-ggplot(d2)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Sexo))+scale_color_manual(values=c("#071e26","#39cccc"))+labs(x="", y="")+theme_bw()ekm3 <-survfit(Surv(time,event)~dukes, data=df2)td1 <- ekm3$time[1:171]td2 <- ekm3$time[172:314]td3 <- ekm3$time[315:length(ekm3$time)]sd1 <- ekm3$surv[1:171]sd2 <- ekm3$surv[172:314]sd3 <- ekm3$surv[315:length(ekm3$surv)]d3 <-data.frame(Tempo =c(td1,td2,td3), `Sobrevivência estimada`=round(c(sd1,sd2,sd3),2), Dukes=c(rep("A-B",length(td1)),rep("C",length(td2)),rep("D", length(td3))))fig3 <-ggplot(d3)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Dukes))+scale_color_manual(values=c("#5639cc","#071e26","#26a4d0"))+labs(x="", y="")+theme_bw()ekm4 <-survfit(Surv(time,event)~charlson, data=df2)tcc1 <- ekm4$time[1:263]tcc2 <- ekm4$time[264:286]tcc3 <- ekm4$time[287:length(ekm4$time)]scc1 <- ekm4$surv[1:263]scc2 <- ekm4$surv[264:286]scc3 <- ekm4$surv[287:length(ekm4$surv)]d4 <-data.frame(Tempo=c(tcc1,tcc2,tcc3), `Sobrevivência estimada`=round(c(scc1,scc2,scc3),2), Charlson =c(rep("0",length(tcc1)),rep("1-2",length(tcc2)),rep("3", length(tcc3))))fig4 <-ggplot(d4)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Charlson))+scale_color_manual(values=c("#5639cc","#071e26","#26a4d0"))+labs(x="", y="")+theme_bw()``````{r}ggarrange(fig1,fig2,fig3,fig4) |>annotate_figure(bottom=text_grob("Tempo até a segunda reinternação (dias)"),left=text_grob("S(t)"))``````{r}valorp <-c(survdiff(Surv(time,event)~chemo, data=df2)$pvalue,survdiff(Surv(time,event)~sex, data=df2)$pvalue,survdiff(Surv(time,event)~dukes, data= df2)$pvalue,survdiff(Surv(time,event)~charlson, data= df2)$pvalue)d <-data.frame(Valorp =round(valorp,3))row.names(d) <-c("chemo", "sex", "dukes", "charlson")kable(d,caption="Valor-p do teste logrank")``````{r}pairwise_survdiff(Surv(time,event)~dukes, data= df2, p.adjust.method ="bonferroni")pairwise_survdiff(Surv(time,event)~charlson, data= df2, p.adjust.method ="bonferroni")```## Tempo até o terceiro eventoConsiderando o tempo até o primeiro evento, a primeira intenação ocorreu no `r min(subset(df3, df3$event == 1)$time)`º dia e a última no `r max(subset(df3, df3$event == 1)$time)`º dia. As curvas de kaplan-Meier podem ser vistas no gráfico a seguir.```{r}ekm1 <-survfit(Surv(time,event)~chemo, data=df3)tc1 <- ekm1$time[1:170]tc2 <- ekm1$time[171:(length(ekm1$time)-1)]sc1 <- ekm1$surv[1:170]sc2 <- ekm1$surv[171:(length(ekm1$time)-1)]d1 <-data.frame(Tempo =c(tc1,tc2), `Sobrevivência estimada`=round(c(sc1,sc2),2), Quimioterapia=c(rep("Não tratado",length(tc1)),rep("Tratado",length(tc2))))fig1 <-ggplot(d1)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Quimioterapia))+scale_color_manual(values=c("#071e26","#26a4d0"))+labs(x="", y="")+theme_bw()ekm2 <-survfit(Surv(time,event)~sex, data=df3)ts1 <- ekm2$time[1:219]ts2 <- ekm2$time[220:(length(ekm2$surv)-1)]ss1 <- ekm2$surv[1:219]ss2 <- ekm2$surv[220:(length(ekm2$surv)-1)]d2 <-data.frame(Tempo =c(ts1,ts2), `Sobrevivência estimada`=round(c(ss1,ss2),2), Sexo=c(rep("Masculino",length(ts1)),rep("Feminino",length(ts2))))fig2 <-ggplot(d2)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Sexo))+scale_color_manual(values=c("#071e26","#39cccc"))+labs(x="", y="")+theme_bw()ekm3 <-survfit(Surv(time,event)~dukes, data=df3)td1 <- ekm3$time[1:170]td2 <- ekm3$time[171:312]td3 <- ekm3$time[313:length(ekm3$time)]sd1 <- ekm3$surv[1:170]sd2 <- ekm3$surv[171:312]sd3 <- ekm3$surv[313:length(ekm3$surv)]d3 <-data.frame(Tempo =c(td1,td2,td3), `Sobrevivência estimada`=round(c(sd1,sd2,sd3),2), Dukes=c(rep("A-B",length(td1)),rep("C",length(td2)),rep("D", length(td3))))fig3 <-ggplot(d3)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Dukes))+scale_color_manual(values=c("#5639cc","#071e26","#26a4d0"))+labs(x="", y="")+theme_bw()ekm4 <-survfit(Surv(time,event)~charlson, data=df3)tcc1 <- ekm4$time[1:266]tcc2 <- ekm4$time[267:285]tcc3 <- ekm4$time[286:length(ekm4$time)]scc1 <- ekm4$surv[1:266]scc2 <- ekm4$surv[267:285]scc3 <- ekm4$surv[286:length(ekm4$surv)]d4 <-data.frame(Tempo=c(tcc1,tcc2,tcc3), `Sobrevivência estimada`=round(c(scc1,scc2,scc3),2), Charlson =c(rep("0",length(tcc1)),rep("1-2",length(tcc2)),rep("3", length(tcc3))))fig4 <-ggplot(d4)+geom_step(aes(x=Tempo, y=`Sobrevivência.estimada`, color=Charlson))+scale_color_manual(values=c("#5639cc","#071e26","#26a4d0"))+labs(x="", y="")+theme_bw()``````{r}ggarrange(fig1,fig2,fig3,fig4) |>annotate_figure(bottom=text_grob("Tempo até a terceira reinternação (dias)"),left=text_grob("S(t)"))``````{r}valorp <-c(survdiff(Surv(time,event)~chemo, data=df3)$pvalue,survdiff(Surv(time,event)~sex, data=df3)$pvalue,survdiff(Surv(time,event)~dukes, data= df3)$pvalue,survdiff(Surv(time,event)~charlson, data= df3)$pvalue)d <-data.frame(Valorp =round(valorp,3))row.names(d) <-c("chemo", "sex", "dukes", "charlson")kable(d,caption="Valor-p do teste logrank")``````{r}pairwise_survdiff(Surv(time,event)~dukes, data= df3, p.adjust.method ="bonferroni")pairwise_survdiff(Surv(time,event)~charlson, data= df3, p.adjust.method ="bonferroni")```# Modelo completo## Abordagem frequentista```{r}sm <-function(op, tipo ="result"){ coef <- op$par se <-SEfromHessian(-1*op$hessian) z <- coef/se p_value <-2*pnorm(abs(z),lower.tail=FALSE) i2.5<- coef+qnorm(0.025)*se i97.5<- coef+qnorm(0.975)*se res <-cbind(coef,exp(coef),exp(-coef),se,z,p_value,i2.5,i97.5)rownames(res) <-c("Intercepto","chemoTreated", "sexFemale", "dukesC", "dukesD", "charlson1-2", "charlson3", "gama","xi")colnames(res) <-c("coef","exp(coef)","exp(-coef)","se(coef)","z","p-value","2.5%","97.5%") result <-as.table(round(res,digits =3))return(result)}``````{r}modelfreq <-coxph(Surv(t.start,t.stop, event)~chemo+sex+dukes+charlson+frailty.gamma(id), data = df)``````{r}X <-cbind(intercepto=rep(1,nrow(df)),model.matrix(modelfreq,data=df)[,-ncol(model.matrix(modelfreq,data=df))])id <- df$idevent <- df$eventtime <- df$timevero <-function(theta){ beta <- theta[1:ncol(X)] gama <- theta[ncol(X)+1] xi <- theta[ncol(X)+2] Dj <-summarise(group_by(df,id),sum(event))[,2] Lambdaj <-summarise(group_by(data.frame(df,d=exp(X%*%beta)*(time^gama)),id),sum(d))[,2] LL <-sum(-(1/xi)*log(xi) -log(gamma(1/xi)) +log(gamma((1/xi)+Dj)) - ((1/xi)+Dj)*log((1/xi)+Lambdaj)) +sum(event*log(gama*(time^(gama-1))*(exp(X%*%beta))))return(LL)}Initial.Values <-c(1,modelfreq$coefficients,0.74,1.3)estfreq <-optim(Initial.Values, vero, hessian =TRUE, method ="L-BFGS-B",control =list(maxit=3000000,fnscale=-1),lower =c(rep(-Inf,ncol(X)),0.01,0.01))estfreq1 <- estfreq```A partir da abordagem frequentista, é possível notar que todas as covariáveis, com exceção da quimioterapia, são significativas ao nível de 5\%. É possível notar ainda que, ao nível de 5\%, há evidências o suficiente para concluir que a variância da fragilidade não é nula e o modelo não se reduz ao modelo para dados independentes.```{r}kable(sm(estfreq))```Tomando o exponencial do seu coeficiente estimado, mantidas fixas as outras covariáveis, tem-se que:- O risco de reinternação dos pacientes do sexo feminino é `r sm(estfreq)[3,2]` vezes o risco dos pacientes do sexo masculino;- O risco de reinternação dos pacientes os pacientes que apresentam estágiotumoral de Dukes A-B é `r sm(estfreq)[4,3]` vezes o risco dos pacientes que apresentam estágio C;- O risco de reinternação dos pacientes os pacientes que apresentam estágiotumoral de Dukes A-B é `r sm(estfreq)[5,3]` vezes o risco dos pacientes que apresentam estágio D;- O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é `r sm(estfreq)[6,3]` vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 1-2;- O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é `r sm(estfreq)[7,3]` vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 3.## Resíduos```{r}#| output: falsetheta <- estfreq1$parbeta <- theta[1:ncol(X)]gama <- theta[ncol(X)+1]xi <- theta[ncol(X)+2]frag <-rgamma(403,1/xi,1/xi)zj <- frag[1]j <-1for (i in2:nrow(df)){if(df$id[i] != df$id[(i-1)]){ j <- j+1 zj <-c(zj,frag[j]) }else{ zj <-c(zj,frag[j]) }}Lambda <-NULLfor (i in1:nrow(df)){ f<-function(t){ lambda <- zj[i]*gama*(t^(gama-1))*(exp(X[i,]%*%beta)) } Lambda <-c(Lambda,integrate(f,0.001,time[i])$value)}resid <-qnorm(Lambda)``````{r}ggplot()+aes(x=seq(1:nrow(df)),y=resid)+geom_point()+labs(x="Índice", y="Resíduos quantílicos")+geom_hline(yintercept =c(3,-3), colour ="blue",linetype ="dashed")+theme_bw() |ggplot() +aes(sample=resid)+stat_qq_point(detrend = T) +stat_qq_band(detrend = T, color='black', fill=NA, size=0.5, linetype="dashed")+geom_hline(yintercept =0, colour ="blue",linetype ="dashed")+geom_vline(xintercept =0, colour ="blue",linetype ="dashed")+labs(x="Quantis teóricos",y="Desvios")+theme_bw()```## Abordagem Bayesiana```{r}X <-cbind(intercepto=rep(1,nrow(df)),model.matrix(modelfreq,data=df)[,-ncol(model.matrix(modelfreq,data=df))])p <-ncol(X)beta0 <-rep(0,p) V0 <-matrix(c(rep(0,p^2)), ncol=p, nrow=p) diag(V0) <-rep(1000,p)c <-0.1Initial.Values <-c(1,modelfreq$coefficients,0.74,1.3)burnin <-1000thin <-200N <- (1000+burnin)*thinparm.names <-as.parm.names(list(beta=rep(0,p),gama=0,xi=0))pos.beta <-grep("beta", parm.names)pos.gama <-grep("gama", parm.names)pos.xi <-grep("xi", parm.names)MyData <-list(J=p, X=X, time=time, id=id, event=event, mon.names="LP", parm.names=parm.names,pos.beta=pos.beta,pos.gama=pos.gama,pos.xi=pos.xi)modelobayes <-function(parm, Data=MyData){ ### Parameters beta <- parm[Data$pos.beta] gama <- parm[Data$pos.gama] xi <- parm[Data$pos.xi]### Log-Prior beta.prior <-dmvn(beta, beta0, V0, log=TRUE) gama.prior <-dexp(gama, c, log=TRUE) xi.prior <-dexp(xi, c, log=TRUE)### Log-Likelihood Dj <-summarise(group_by(df,id),sum(event))[,2] Lambdaj <-summarise(group_by(data.frame(df,d=exp(X%*%beta)*(time^gama)),id),sum(d))[,2] LL <-sum(-(1/xi)*log(xi) -log(gamma(1/xi)) +log(gamma((1/xi)+Dj)) - ((1/xi)+Dj)*log((1/xi)+Lambdaj)) +sum(event*log(gama*(time^(gama-1))*(exp(X%*%beta))))### Log-Posterior LP <- LL + beta.prior + gama.prior + xi.prior Modelout <-list(LP=LP, Dev=-2*LL, Monitor=LP, yhat =exp(X*beta), parm=parm)return(Modelout)}``````{r}#| output: false#modelbayes <- LaplacesDemon(modelobayes, Data=MyData, Initial.Values, Covar=NULL, #Iterations=N, Status=N/5, Thinning=thin, Algorithm="AMWG", Specs=NULL)load("G:/Meu Drive/TCC/Scripts/model_bayes_weibull1.RData")modelbayes1 <- modelbayes``````{r}beta1 <-apply(modelbayes$Posterior1,2,mlv)vero1 <-vero(beta1)prior1 <-dmvn(beta1[1:ncol(X)], beta0, V0, log=TRUE)+dexp(beta1[(ncol(X)+1)], c, log=TRUE)+dexp(beta1[(ncol(X)+2)], c, log=TRUE)```Utilizando o algoritmo Adaptative Metropolis-within-Gibbs, obteve-se uma taxa de aceitação de `r round(modelbayes$Acceptance.Rate*100,2)`\%, próxima do considerado ótimo.```{r}dic1 <-data.frame(all=round(modelbayes$DIC1,2), stationary=round(modelbayes$DIC1,2))colnames(dic1) <-c("Todas as amostras","Estacionárias")rownames(dic1) <-c("Dbar", "pd", "DIC")``````{r}mode <-cbind(c(round(apply(modelbayes$Posterior1,2,mlv),3),Deviance=round(mlv(modelbayes$Deviance),3),LP="-"))d <-data.frame(mode = mode,'exp(coef)'=c(round(exp(as.numeric(mode[1:9])),3), "-","-"),'exp(-coef)'=c(round(exp(-as.numeric(mode[1:9])),3),"-","-"),round(modelbayes$Summary1,3))```Tomando o exponencial do seu coeficiente estimado, mantidas fixas as outras covariáveis, tem-se que:- O risco de reinternação dos pacientes do sexo feminino é `r d[3,2]` vezes o risco dos pacientes do sexo masculino;- O risco de reinternação dos pacientes os pacientes que apresentam estágiotumoral de Dukes A-B é `r d[4,3]` vezes o risco dos pacientes que apresentam estágio C;- O risco de reinternação dos pacientes os pacientes que apresentam estágiotumoral de Dukes A-B é `r d[5,3]` vezes o risco dos pacientes que apresentam estágio D;- O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é `r d[6,3]` vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 1-2;- O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é `r d[7,3]` vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 3;```{r}kable(d, caption ="Tabela: Medidas resumo das distribuições dos parâmetros da distribuição a posteriori, do desvio e do LP.")```A partir dos gráficos a seguir e da Tabela, é possível notar que as distribuições a posteriori dos betas, gama e xi são simétricas. É possível notar também há uma notável estacionariedade, evidenciada pela ausência de tendências nos gráficos. O que indica que as amostras estão sendo colhidas de forma aleatória da mesma distribuição alvo. Os correlogramas exibidos revelam uma escassa presença de autocorrelações significativas, sugerindo que as amostras podem ser consideradas independentes. Ou seja, as amostras geradas não são fortemente influenciadas pelas amostras geradas anteriormente. Os desvios sugerem que a maioria das verossimilhanças foram maximizadas durante o processo, já que a distribuição é assimétrica à direita. Além disso, o correlograma não apresenta indícios de autocorrelação nos desvios das amostras.Os logaritmos da distribuição a posteriori, das amostras indicam a sua boa adequação, uma vez que a sua distribuição é assimétrica à esquerda e o correlograma apresenta apresenta apenas três correlações significativas. ::: panel-tabset### Beta 0```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,1])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,1]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[0]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,1])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,1],60)+labs(title="")+theme_bw()```### Beta 1```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,2])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,2]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[1]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,2])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,2],60)+labs(title="")+theme_bw()```### Beta 2```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,3])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,3]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[2]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,3])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,3],60)+labs(title="")+theme_bw()```### Beta 3```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,4])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,4]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[3]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,4])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,4],60)+labs(title="")+theme_bw()```### Beta 4```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,5])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,5]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[4]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,5])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,5],60)+labs(title="")+theme_bw()```### Beta 5```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,6])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,6]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[5]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,6])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,6],60)+labs(title="")+theme_bw()```### Beta 6```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,7])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,7]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[6]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,7])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,7],60)+labs(title="")+theme_bw()```### Gama```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,8])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,8]), color ="#50a3fc")+labs(x="Iteração", y=expression(gamma))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,8])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,8],60)+labs(title="")+theme_bw()```### xi```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,9])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,9]), color ="#50a3fc")+labs(x="Iteração", y=expression(xi))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,9])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,9],60)+labs(title="")+theme_bw()```### Deviance```{r}ggplot()+aes(x=1:length(modelbayes$Deviance), y=modelbayes$Deviance)+geom_line() +geom_hline(yintercept=mean(modelbayes$Deviance), color ="#50a3fc")+labs(x="Iteração", y="Desvio")+theme_bw()+ggplot()+aes(x=modelbayes$Deviance)+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Deviance,60)+labs(title="")+theme_bw()```### LP```{r}ggplot()+aes(x=1:length(modelbayes$Monitor), y=modelbayes$Monitor)+geom_line() +geom_hline(yintercept=mean(modelbayes$Monitor), color ="#50a3fc")+labs(x="Iteração", y="LP")+theme_bw()+ggplot()+aes(x=modelbayes$Monitor)+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Monitor,60)+labs(title="")+theme_bw()```:::## Resíduos```{r}#| output: falsetheta <-apply(modelbayes$Posterior1,2,mlv)beta <- theta[1:ncol(X)]gama <- theta[ncol(X)+1]xi <- theta[ncol(X)+2]frag <-rgamma(403,1/xi,1/xi)zj <- frag[1]j <-1for (i in2:nrow(df)){if(df$id[i] != df$id[(i-1)]){ j <- j+1 zj <-c(zj,frag[j]) }else{ zj <-c(zj,frag[j]) }}Lambda <-NULLfor (i in1:nrow(df)){ f<-function(t){ lambda <- zj[i]*gama*(t^(gama-1))*(exp(X[i,]%*%beta)) } Lambda <-c(Lambda,integrate(f,0.001,time[i])$value)}resid <-qnorm(Lambda)``````{r}ggplot()+aes(x=seq(1:nrow(df)),y=resid)+geom_point()+labs(x="Índice", y="Resíduos quantílicos")+geom_hline(yintercept =c(3,-3), colour ="blue",linetype ="dashed")+theme_bw() |ggplot() +aes(sample=resid)+stat_qq_point(detrend = T) +stat_qq_band(detrend = T, color='black', fill=NA, size=0.5, linetype="dashed")+geom_hline(yintercept =0, colour ="blue",linetype ="dashed")+geom_vline(xintercept =0, colour ="blue",linetype ="dashed")+labs(x="Quantis teóricos",y="Desvios")+theme_bw()```# Modelo reduzido## Abordagem frequentista```{r}df <-read.csv("readmission.csv")|>mutate(chemo =NULL,sex =relevel(factor(sex), "Male"),dukes =relevel(factor(dukes), "A-B"),charlson =relevel(factor(charlson), "0"),death =NULL,time = t.stop) ``````{r}sm <-function(op, tipo ="result"){ coef <- op$par se <-SEfromHessian(-1*op$hessian) z <- coef/se p_value <-2*pnorm(abs(z),lower.tail=FALSE) i2.5<- coef+qnorm(0.025)*se i97.5<- coef+qnorm(0.975)*se res <-cbind(coef,exp(coef),exp(-coef),se,z,p_value,i2.5,i97.5)rownames(res) <-c("Intercepto", "sexFemale", "dukesC", "dukesD", "charlson1-2", "charlson3", "gama","xi")colnames(res) <-c("coef","exp(coef)","exp(-coef)","se(coef)","z","p-value","2.5%","97.5%") result <-as.table(round(res,digits =3))return(result)}``````{r}modelfreq <-coxph(Surv(t.start,t.stop, event)~sex+dukes+charlson+frailty.gamma(id), data = df)``````{r}X <-cbind(intercepto=rep(1,nrow(df)),model.matrix(modelfreq,data=df)[,-ncol(model.matrix(modelfreq,data=df))])id <- df$idevent <- df$eventtime <- df$timevero <-function(theta){ beta <- theta[1:ncol(X)] gama <- theta[ncol(X)+1] xi <- theta[ncol(X)+2] Dj <-summarise(group_by(df,id),sum(event))[,2] Lambdaj <-summarise(group_by(data.frame(df,d=exp(X%*%beta)*(time^gama)),id),sum(d))[,2] LL <-sum(-(1/xi)*log(xi) -log(gamma(1/xi)) +log(gamma((1/xi)+Dj)) - ((1/xi)+Dj)*log((1/xi)+Lambdaj)) +sum(event*log(gama*(time^(gama-1))*(exp(X%*%beta))))return(LL)}Initial.Values <-c(1,modelfreq$coefficients,0.74,1.3)estfreq <-optim(Initial.Values, vero, hessian =TRUE, method ="L-BFGS-B",control =list(maxit=3000000,fnscale=-1),lower =c(rep(-Inf,ncol(X)),0.01,0.01))estfreq2 <- estfreq```A partir da abordagem frequentista, é possível notar que todas as covariáveis são significativas ao nível de 5\%. É possível notar ainda que, ao nível de 5\%, há evidências o suficiente para concluir que a variância da fragilidade não é nula e o modelo não se reduz ao modelo para dados independentes.```{r}kable(sm(estfreq))```Tomando o exponencial do seu coeficiente estimado, mantidas fixas as outras covariáveis, tem-se que:- O risco de reinternação dos pacientes do sexo feminino é `r sm(estfreq)[2,2]` vezes o risco dos pacientes do sexo masculino;- O risco de reinternação dos pacientes os pacientes que apresentam estágiotumoral de Dukes A-B é `r sm(estfreq)[3,3]` vezes o risco dos pacientes que apresentam estágio C;- O risco de reinternação dos pacientes os pacientes que apresentam estágiotumoral de Dukes A-B é `r sm(estfreq)[4,3]` vezes o risco dos pacientes que apresentam estágio D;- O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é `r sm(estfreq)[5,3]` vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 1-2;- O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é `r sm(estfreq)[6,3]` vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 3.## Resíduos```{r}#| output: falsetheta <- estfreq2$parbeta <- theta[1:ncol(X)]gama <- theta[ncol(X)+1]xi <- theta[ncol(X)+2]frag <-rgamma(403,1/xi,1/xi)zj <- frag[1]j <-1for (i in2:nrow(df)){if(df$id[i] != df$id[(i-1)]){ j <- j+1 zj <-c(zj,frag[j]) }else{ zj <-c(zj,frag[j]) }}Lambda <-NULLfor (i in1:nrow(df)){ f<-function(t){ lambda <- zj[i]*gama*(t^(gama-1))*(exp(X[i,]%*%beta)) } Lambda <-c(Lambda,integrate(f,0.001,time[i])$value)}resid <-qnorm(Lambda)``````{r}ggplot()+aes(x=seq(1:nrow(df)),y=resid)+geom_point()+labs(x="Índice", y="Resíduos quantílicos")+#geom_hline(yintercept = c(3,-3), colour = "blue",linetype = "dashed")+theme_bw() |ggplot() +aes(sample=resid)+stat_qq_point(detrend = T) +stat_qq_band(detrend = T, color='black', fill=NA, size=0.5, linetype="dashed")+geom_hline(yintercept =0, colour ="blue",linetype ="dashed")+geom_vline(xintercept =0, colour ="blue",linetype ="dashed")+labs(x="Quantis teóricos",y="Desvios")+theme_bw()```## Abordagem Bayesiana```{r}X <-cbind(intercepto=rep(1,nrow(df)),model.matrix(modelfreq,data=df)[,-ncol(model.matrix(modelfreq,data=df))])p <-ncol(X)beta0 <-rep(0,p) V0 <-matrix(c(rep(0,p^2)), ncol=p, nrow=p) diag(V0) <-rep(1000,p)c <-0.1Initial.Values <-c(1,modelfreq$coefficients,0.74,1.3)burnin <-1000thin <-200N <- (1000+burnin)*thinparm.names <-as.parm.names(list(beta=rep(0,p),gama=0,xi=0))pos.beta <-grep("beta", parm.names)pos.gama <-grep("gama", parm.names)pos.xi <-grep("xi", parm.names)MyData <-list(J=p, X=X, time=time, id=id, event=event, mon.names="LP", parm.names=parm.names,pos.beta=pos.beta,pos.gama=pos.gama,pos.xi=pos.xi)modelobayes <-function(parm, Data=MyData){ ### Parameters beta <- parm[Data$pos.beta] gama <- parm[Data$pos.gama] xi <- parm[Data$pos.xi]### Log-Prior beta.prior <-dmvn(beta, beta0, V0, log=TRUE) gama.prior <-dexp(gama, c, log=TRUE) xi.prior <-dexp(xi, c, log=TRUE)### Log-Likelihood Dj <-summarise(group_by(df,id),sum(event))[,2] Lambdaj <-summarise(group_by(data.frame(df,d=exp(X%*%beta)*(time^gama)),id),sum(d))[,2] LL <-sum(-(1/xi)*log(xi) -log(gamma(1/xi)) +log(gamma((1/xi)+Dj)) - ((1/xi)+Dj)*log((1/xi)+Lambdaj)) +sum(event*log(gama*(time^(gama-1))*(exp(X%*%beta))))### Log-Posterior LP <- LL + beta.prior + gama.prior + xi.prior Modelout <-list(LP=LP, Dev=-2*LL, Monitor=LP, yhat =exp(X*beta), parm=parm)return(Modelout)}``````{r}#| output: false#modelbayes <- LaplacesDemon(modelobayes, Data=MyData, Initial.Values, Covar=NULL, #Iterations=N, Status=N/5, Thinning=thin, Algorithm="AMWG", Specs=NULL)load("G:/Meu Drive/TCC/Scripts/model_bayes_weibull2.RData")modelbayes2 <- modelbayes``````{r}beta2 <-apply(modelbayes$Posterior1,2,mlv)vero2 <-vero(beta2)prior2 <-dmvn(beta2[1:ncol(X)], beta0, V0, log=TRUE)+dexp(beta2[(ncol(X)+1)], c, log=TRUE)+dexp(beta2[(ncol(X)+2)], c, log=TRUE)```Utilizando o algoritmo Adaptative Metropolis-within-Gibbs, obteve-se uma taxa de aceitação de `r round(modelbayes$Acceptance.Rate*100,2)`\%, próxima do considerado ótimo.```{r}dic2 <-data.frame(all=round(modelbayes$DIC1,2), stationary=round(modelbayes$DIC1,2))colnames(dic2) <-c("Todas as amostras","Estacionárias")rownames(dic2) <-c("Dbar", "pd", "DIC")``````{r}mode <-cbind(c(round(apply(modelbayes$Posterior1,2,mlv),3),Deviance=round(mlv(modelbayes$Deviance),3),LP="-"))d <-data.frame(mode = mode,'exp(coef)'=c(round(exp(as.numeric(mode[1:8])),3), "-","-"),'exp(-coef)'=c(round(exp(-as.numeric(mode[1:8])),3),"-","-"),round(modelbayes$Summary1,3))```Tomando o exponencial do seu coeficiente estimado, mantidas fixas as outras covariáveis, tem-se que:- O risco de reinternação dos pacientes do sexo feminino é `r d[2,2]` vezes o risco dos pacientes do sexo masculino;- O risco de reinternação dos pacientes os pacientes que apresentam estágiotumoral de Dukes A-B é `r d[3,3]` vezes o risco dos pacientes que apresentam estágio C;- O risco de reinternação dos pacientes os pacientes que apresentam estágiotumoral de Dukes A-B é `r d[4,3]` vezes o risco dos pacientes que apresentam estágio D;- O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é `r d[5,3]` vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 1-2;- O risco de reinternação dos pacientes os pacientes que apresentam índice de comorbidade de Charlson 0 é `r d[6,3]` vezes o risco dos pacientes os pacientes que apresentam índice de comorbidade 3;```{r}kable(d, caption ="Tabela: Medidas resumo das distribuições dos parâmetros da distribuição a posteriori, do desvio e do LP.")```A partir dos gráficos a seguir e da Tabela, é possível notar que as distribuições a posteriori dos betas, gama e xi são simétricas. É possível notar também há uma notável estacionariedade, evidenciada pela ausência de tendências nos gráficos. O que indica que as amostras estão sendo colhidas de forma aleatória da mesma distribuição alvo. Os correlogramas exibidos revelam uma escassa presença de autocorrelações significativas, sugerindo que as amostras podem ser consideradas independentes. Ou seja, as amostras geradas não são fortemente influenciadas pelas amostras geradas anteriormente. Os desvios sugerem que a maioria das verossimilhanças foram maximizadas durante o processo, já que a distribuição é assimétrica à direita. Além disso, o correlograma não apresenta indícios de autocorrelação nos desvios das amostras.Os logaritmos da distribuição a posteriori, das amostras indicam a sua boa adequação, uma vez que a sua distribuição é assimétrica à esquerda e o correlograma apresenta apresenta apenas três correlações significativas. ::: panel-tabset### Beta 0```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,1])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,1]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[0]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,1])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,1],60)+labs(title="")+theme_bw()```### Beta 1```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,2])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,2]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[1]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,2])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,2],60)+labs(title="")+theme_bw()```### Beta 2```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,3])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,3]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[2]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,3])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,3],60)+labs(title="")+theme_bw()```### Beta 3```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,4])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,4]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[3]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,4])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,4],60)+labs(title="")+theme_bw()```### Beta 4```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,5])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,5]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[4]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,5])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,5],60)+labs(title="")+theme_bw()```### Beta 5```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,6])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,6]), color ="#50a3fc")+labs(x="Iteração", y=expression(beta[5]))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,6])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,6],60)+labs(title="")+theme_bw()```### Gama```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,7])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,7]), color ="#50a3fc")+labs(x="Iteração", y=expression(gamma))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,7])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,7],60)+labs(title="")+theme_bw()```### xi```{r}ggplot()+aes(x=1:nrow(modelbayes$Posterior1), y=modelbayes$Posterior1[,8])+geom_line() +geom_hline(yintercept=mean(modelbayes$Posterior1[,8]), color ="#50a3fc")+labs(x="Iteração", y=expression(xi))+theme_bw()+ggplot()+aes(x=modelbayes$Posterior1[,8])+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Posterior1[,8],60)+labs(title="")+theme_bw()```### Deviance```{r}ggplot()+aes(x=1:length(modelbayes$Deviance), y=modelbayes$Deviance)+geom_line() +geom_hline(yintercept=mean(modelbayes$Deviance), color ="#50a3fc")+labs(x="Iteração", y="Desvio")+theme_bw()+ggplot()+aes(x=modelbayes$Deviance)+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Deviance,60)+labs(title="")+theme_bw()```### LP```{r}ggplot()+aes(x=1:length(modelbayes$Monitor), y=modelbayes$Monitor)+geom_line() +geom_hline(yintercept=mean(modelbayes$Monitor), color ="#50a3fc")+labs(x="Iteração", y="LP")+theme_bw()+ggplot()+aes(x=modelbayes$Monitor)+geom_density(colour="#003e62", fill="#50a3fc")+labs(x="Densidade", y="Valores")+theme_bw() +ggAcf(modelbayes$Monitor,60)+labs(title="")+theme_bw()```:::## Resíduos```{r}#| output: falsetheta <-apply(modelbayes$Posterior1,2,mlv)beta <- theta[1:ncol(X)]gama <- theta[ncol(X)+1]xi <- theta[ncol(X)+2]frag <-rgamma(403,1/xi,1/xi)zj <- frag[1]j <-1for (i in2:nrow(df)){if(df$id[i] != df$id[(i-1)]){ j <- j+1 zj <-c(zj,frag[j]) }else{ zj <-c(zj,frag[j]) }}Lambda <-NULLfor (i in1:nrow(df)){ f<-function(t){ lambda <- zj[i]*gama*(t^(gama-1))*(exp(X[i,]%*%beta)) } Lambda <-c(Lambda,integrate(f,0.001,time[i])$value)}resid <-qnorm(Lambda)``````{r}ggplot()+aes(x=seq(1:nrow(df)),y=resid)+geom_point()+labs(x="Índice", y="Resíduos quantílicos")+#geom_hline(yintercept = c(3,-3), colour = "blue",linetype = "dashed")+theme_bw() |ggplot() +aes(sample=resid)+stat_qq_point(detrend = T) +stat_qq_band(detrend = T, color='black', fill=NA, size=0.5, linetype="dashed")+geom_hline(yintercept =0, colour ="blue",linetype ="dashed")+geom_vline(xintercept =0, colour ="blue",linetype ="dashed")+labs(x="Quantis teóricos",y="Desvios")+theme_bw()``````{r}#dens <- function(theta){# beta <- theta[1:ncol(X)]# gama <- theta[ncol(X)+1]# xi <- theta[ncol(X)+2]# # Dj <- summarise(group_by(df,id),sum(event))[,2]# Lambdaj <- summarise(group_by(data.frame(df,d=exp(X%*%beta)*(time^gama)),id),sum(d))[,2]# dj <- summarise(group_by(data.frame(df,d=event*log(gama*(time^(gama-1))*(exp(X%*%beta)))),id),sum(d))[,2]# LF <- -(1/xi)*log(xi) -log(gamma(1/xi)) + log(gamma((1/xi)+Dj)) - ((1/xi)+Dj)*log((1/xi)+Lambdaj) + dj# return(lF)#}```# Comparação de modelos```{r}RV <--2*(-estfreq1$value+estfreq2$value)pRV <-pchisq(RV,1,lower.tail=FALSE)```Na abordagem frequentista, o teste da razão de verossimilhanças resultou num valor-p de `r pRV`. Dessa forma, não rejeitamos a hipótese nula ao nível de 5\%. Ou seja, há evidências o suficiente para concluir que o modelo reduzido é mais adequado que o modelo completo. Além disso, o AIC do modelo reduzido é menor, indicando a sua melhor adequação.```{r}AIC_1 <--2*estfreq1$value+2*7BIC_1 <--2*estfreq1$value+5*log(nrow(df))AIC_2 <--2*estfreq2$value+2*6BIC_2 <--2*estfreq2$value+5*log(nrow(df))d <-data.frame(c(AIC_1,BIC_1),c(AIC_2,BIC_2))colnames(d) <-c("Modelo completo", "Modelo reduzido")row.names(d) <-c("AIC","BIC")kable(d, caption="Tabela: AIC e BIC dos modelos completo e reduzido para os dados de cancer colorretal.")```Na abordagem Bayesiana, o DIC para o modelo reduzido é menor que o do modelo completo, indicando a sua melhor adequação.```{r}d <-data.frame(dic1[,1],dic2[,1])rownames(d) <-c("Dbar", "pd", "DIC")colnames(d) <-c("Modelo completo","Modelo reduzido")kable(d, caption="Tabela: DIC das amostras a posteriori dos modelos completo e reduzido para os dados de cancer colorretal.")``````{r}k <- (vero1*prior2)/(vero2*prior1)```O fator de Bayes é `r round(k,2)`, isso indica que há evidências a favor do modelo reduzido. Além disso, o modelo com todas as covariáveis precisou de `r round(modelbayes1$Minutes/60,2)`h para gerar todas as amostras, enquanto o modelo sem a covariável quimioterapia precisou de `r round(modelbayes2$Minutes/60,2)`h.